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Abstract 

Early  identification  of  impending  illness  during  widespread  exposure  to  a  pathogenic  agent 
offers  a  potential  means  to  initiate  treatment  during  a  timeframe  when  it  would  be  most  likely 
to  be  effective  and  has  the  potential  to  identify  novel  therapeutic  strategies.  The  latter  could 
be  critical,  especially  as  antibiotic  resistance  is  becoming  widespread.  In  order  to  examine 
pre-symptomatic  illness,  African  green  monkeys  were  challenged  intranasally  with  aerosol¬ 
ized  Yersinia  pestis  strain  C092  and  blood  samples  were  collected  in  short  intervals  from 
45  m  till  42  h  post-exposure.  Presenting  one  of  the  first  genomic  investigations  of  a  NHP 
model  challenged  by  pneumonic  plague,  whole  genome  analysis  was  annotated  in  silico 
and  validated  by  qPCR  assay.  Transcriptomic  profiles  of  blood  showed  early  perturbation 
with  the  number  of  differentially  expressed  genes  increasing  until  24  h.  By  then,  Y.  pestis 
had  paralyzed  the  host  defense,  as  suggested  by  the  functional  analyses.  Early  activation 
of  the  apoptotic  networks  possibly  facilitated  the  pathogen  to  overwhelm  the  defense  mech¬ 
anisms,  despite  the  activation  of  the  pro-inflammatory  mechanism,  toll-like  receptors  and 
microtubules  at  the  port-of-entry.  The  overexpressed  transcripts  encoding  an  early  pro- 
inflammatory  response  particularly  manifested  in  active  lymphocytes  and  ubiquitin  networks 
were  a  potential  deviation  from  the  rodent  models,  which  needs  further  verification.  In  sum¬ 
mary,  the  present  study  recognized  a  pattern  of  Y.  pestis  pathogenesis  potentially  more 
applicable  to  the  human  system.  Independent  validation  using  the  complementary  omics 
approach  with  comprehensive  evaluation  of  the  organs,  such  as  lungs  which  showed  early 
bacterial  infection,  is  essential. 
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Introduction 

Extracellular  proliferation  and  rapid  propagation  of  Yersinia  distinguishes  its  pathogenesis 
from  other  enteropathogenic  agents,  such  as  Salmonella  and  Shigella  [1-3].  The  capability  of 
Yersinia  pestis  (Y.  pestis)  to  rapidly  multiply,  coupled  with  its  atypical  symptoms  at  the  initial 
stage  of  infection,  presents  daunting  challenges  for  initiating  immediate  therapeutic  action  and 
appropriate  preventative  care.  Y.  pestis,  the  causative  agent  of  a  pneumonic  plague,  manifests 
early  flu-like  symptoms  that  evolve  into  chest  pain,  shortness  of  breath,  cough,  and  bloody  or 
watery  sputum.  If  untreated  within  the  first  20  hours  of  onset  of  symptoms,  it  is  almost  always 
fatal  due  to  respiratory  failure  and/or  circulatory  collapse,  coagulopathy,  and  hemorrhage  [4], 
Its  alarmingly  high  mortality  rate,  even  under  antimicrobial  treatment,  and  its  easy  transmis¬ 
sion  and  aerosol  characteristics  continually  consider  Y.  pestis  as  a  major  threat  to  public  health 
[5]  and  potential  agent  of  bioweapon  [6].  Recent  fatality  caused  by  Y.  pestis  infection  has  been 
attributed  to  the  clinical  misdiagnosis  and  misidentification  by  multiple  automated  platforms 
further  underlines  the  need  for  a  reliable  early  detection  technique  [7]. 

Y.  pestis  is  a  rod-shaped  Gram-negative  intracellular  bacterium  that,  once  inhaled,  first 
encounters  phagocytosis  mediated  by  macrophages  and  neutrophils  [2,8].  Surviving  pathogens 
inside  the  macrophages  proliferate,  continually  expressing  various  virulence  factors  such  as 
Yops,  the  FI  antigen,  and  the  V  antigen.  Once  released  into  the  blood  stream,  the  extracellular 
pathogens  can  rapidly  gain  access  to  many  organs  and  systematically  modulate  the  innate  and 
adaptive  immune  responses.  Pathological  analysis  reported  a  rapid  dissemination  of  Y.  pestis  in 
murine  lungs  within  1  h  after  the  intranasal  challenge  that  developed  into  a  splenic  infection 
within  36  h  [9].  An  ex  vivo  bioluminescence  study  found  Y.  pestis  proliferated  in  lungs  1  d  after 
intranasal  challenge;  within  24  h,  pathogens  had  spread  to  the  spleen  and  liver,  and  by  72  h  the 
entire  body  was  infected  [10],  Interestingly,  no  viable  pathogens  were  detected  in  rodent  [11] 
and  primate  [12]  blood  until  72  h  post-exposure. 

A  comprehensive  investigation  of  the  early  host-pathogen  dynamics  would  be  of  a  great 
importance,  particularly  in  regard  to  biological  warfare  and  clinical  profiles.  The  study  of  the 
infection  trajectory  focused  on  the  important  clinical  and  molecular  landmarks  of  disease  pro¬ 
gression  could  be  highly  insightful,  yet  is  largely  lacking.  The  in  vivo  chain  of  events  following 
the  inhalation  of  Y.  pestis  is  still  poorly  understood,  especially  in  the  case  of  large  mammals, 
thereby  presenting  a  critical  knowledge  gap  [13,14]. 

In  the  past,  omics  approaches  were  used  to  interrogate  the  longitudinal  murine  pneumonic 
plague  model  [15,16],  while  increasing  interest  has  been  shown  on  its  early  pathogenesis  [17- 
19].  Similar  to  the  distinguishing  signature  of  the  bubonic  plague  model  [20,21],  pneumonic 
plague  in  rodents  manifested  a  biphasic  mode  of  temporal  progression  in  lungs,  kidney  and 
spleen  tissues  [9,22].  Here,  the  early  pre-inflammatory  phase  was  found  to  have  reduced  tran- 
scriptomic  recruitment  until  12  h  post  infection,  followed  by  a  pro-inflammatory  augmenta¬ 
tion  that  resulted  in  a  larger  number  of  genes  dysregulated  after  48  h  post-challenge  [15].  The 
majority  of  the  inflammation  markers  was  overexpressed  in  blood,  lungs,  spleen,  liver  and 
heart  at  72  h  post- infection  [1 1,15].  A  similar  observation  was  also  reported  in  rats  challenged 
intranasally  with  Y.  pestis  [23]. 

Investigation  of  the  nonhuman  primates  that  share  phylogenetic  proximity  to  humans 
becomes  essential  to  bridge  the  lack  of  clinical  knowledge  about  Y.  pestis.  In  the  present  study, 
we  characterized  the  longitudinal  progression  of  pneumonic  plague  in  African  green  monkeys 
( Chlorocebus  aethiops )  exposed  to  Y.  pestis  strain  C092  via  inhalation  challenge.  Blood  draws 
were  designed  with  shorter  intervals,  starting  at  45  m  post-exposure,  to  take  a  more  detailed 
look  at  the  disease  trajectory.  Layton  et  al.,  [12]  characterized  this  animal  model  as  a  viable  sub¬ 
stitute  for  human  testing  according  to  the  FDA  ‘Animal  Rule’.  They  observed  that  the  bacterial 
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load  became  detectable  in  blood  in  a  little  less  than  three  days.  On  average,  the  animals  sur¬ 
vived  another  day  or  two.  By  then,  the  pathogen  was  widely  disseminated  in  tissues  such  as  the 
lymph  nodes,  lungs,  liver,  and  spleen  [12],  Following  the  published  animal  handling  guidelines 
[12,24],  we  collected  blood  at  eight  time  points  from  the  moribund  animals,  starting  from  45  m 
until  42  h.  High  throughput  transcriptomic  outcome  was  annotated  in  silico  and  validated  by 
qPCR  assay. 

Materials  and  Methods 

Ethics  Statement 

All  animal  experiments  were  approved  by  the  Institutional  Animal  Care  and  Use  Committee 
(IACUC)  at  the  Walter  Reed  Army  Institute  of  Research  (WRAIR),  Silver  Spring,  MD,  and 
were  performed  in  a  facility  accredited  by  the  Association  for  the  Assessment  and  Accredita¬ 
tion  of  Laboratory  Animal  Care  International  (AAALAC).  The  approved  protocol  was  POOl- 
08:  Systems  biology  studies  to  identify  host  indicators  /therapeutic  targets  at  early  time  periods 
post  exposure  to  Y.pestis  (C092)  in  African  Green  Monkeys  ( Chlorocebus  aethiops) 

Animals 

The  animals  were  obtained  from  the  primate  colony  of  WRAIR.  Adult  male  African  green 
monkeys  ( Chlorocebus  aethiops)  of  4. 8-7.0  kg  body  weight  were  TB-negative,  SIV-negative, 
and  SRV-negative  with  serum  titers  of  anti-Yersinia  IgG  antibodies  <1:160.  Only  monkeys  in 
good  health  as  shown  by  physical  examination,  complete  blood  counts,  and  serum  chemistry 
profiles  were  used.  Monkeys  were  individually  housed  in  3  x  2.6  x  2-ft.  squeeze  cages  in  free¬ 
standing  enclosures  equipped  with  standard  enrichments,  exposed  to  ambient  environmental 
conditions  inside  an  Animal  Biosafety  Level  3  (ABSL-3)  containment  laboratory,  which 
approximated  a  12: 12-hour  light/dark  cycle  with  temperatures  between  25°C  and  30°C.  Mon¬ 
keys  were  fed  a  standard  laboratory  primate  chow  (TekLad,  Madison,  WI)  with  water  provided 
ad  libitum  by  Lixit  valve  (Lixit  Corp.,  Napa,  CA).  Veterinary  care  was  provided  twice  a  day,  but 
not  more  than  12  h  interval  to  the  healthy  animals.  Once  exposed,  the  interval  time  was 
reduced  to  not  longer  than  2  h. 

Overall  experimental  design  (SI  Fig) 

Following  a  published  report  [12],  the  present  study  was  designed  to  expose  the  nonhuman  pri¬ 
mates  with  a  target  dose  of  100  ±  50  LD50  aerosolized  Y.pestis  strain  C092  under  ABSL-3  con¬ 
ditions.  Seventy-two  hours  prior  to  aerosol  exposure,  all  animals  were  moved  into  the  ABSL-3 
for  acclimatization.  As  described  in  SI  Fig,  initial  blood  draws  were  carried  out  24  h  before  Y. 
pestis  aerosol  challenge.  The  animals  fasted  for  last  6  h  and  were  anesthetized  using  4  mg/kg 
Telazol  (Fort  Dodge  Animal  Health,  Fort  Dodge,  IA);  after  15  m,  they  were  aerosol  challenged. 
Post- exposure  clinical  illness  and  disseminated  infection  were  evaluated  by  direct  clinical 
observations.  All  animals  within  each  group  were  exposed  on  the  same  day  at  30-minute  inter¬ 
vals,  and  subsequent  blood  draws  were  conducted  at  45  m,  6  h,  9  h,  12  h,  18  h,  24  h,  32  h  and 
42  h  post-exposure. 

Except  for  the  45  m  blood  draw,  the  remaining  blood  draws  were  terminal  with  the  animals 
immediately  euthanized  to  collect  organs.  For  the  45  m  and  6  h  blood  draws,  we  used  the  same 
cohort.  The  blood  draws  were  carried  out  at  the  same  time  on  the  scheduled  days  to  minimize 
molecular  variances  attributed  to  the  circadian  cycle.  All  animals  were  humanely  euthanized 
by  first  inducing  a  deep  plane  of  anesthesia  with  intramuscular  administration  of  9  mg/kg  Tela¬ 
zol.  Euthasol  (390  mg/ml  pentobarbital  and  50  mg/ml  phenytoin  (Virbac  Corp.,  Fort  Worth, 
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TX,  IC))  was  administered  at  a  dose  of  60  mg/kg  pentobarbital  sodium  using  an  18-gauge  nee¬ 
dle.  The  animal  was  auscultated  to  determine  the  absence  of  heartbeat  and  breathing.  Lethality 
attributed  to  the  given  pathogenic  load  was  tested  in  two  animals  which  met  their  end  points 
around  78  h  post-exposure. 

The  IACUC  approved  PO01-08  protocol  listed  nine  clinical  parameters  (SI  Table);  any  ani¬ 
mal  that  met  3  or  more  of  those  clinical  parameters  or  was  moribund,  cachectic,  or  unable  to 
obtain  food  or  water  would  be  humanely  euthanized.  Based  on  this  guideline,  none  of  the  two 
marked  for  lethality  test  did  not  meet  sufficient  clinical  parameters  of  euthanasia  during  the 
last  health  checkup.  However,  within  just  2  hours  between  two  successive  animal  health  inspec¬ 
tions  one  animal  succumbed  to  the  disease;  unfortunately  our  very  competent  team  did  not  see 
the  signs  of  imminent  mortality.  The  other  one  marked  for  the  lethality  test  was  euthanized  at 
this  time  after  consulting  with  staff  Veterinarian.  To  note:  the  NHPs  are  notoriously  able  to 
hide  their  illness.  The  pathogenesis  of  pneumonic  plague  in  this  animal  was  very  rapid  and 
even  though  the  veterinary  staff  took  maximum  care  making  observations,  the  one  animal  suc¬ 
cumbed  quickly  between  the  two  observation  time  points  (just  2  h).  This  is  true  in  humans  as 
well  and  that  is  why  the  FDA  named  the  African  Green  NHP  as  the  “one  animal  rule”  to  repli¬ 
cate  the  progression  of  illness  upon  infection  with  virulent  Y.  pestis. 

Bacterial  inoculation  and  aerosol  delivery 

Cultures  of  Y.  pestis  C092  strain  were  grown  overnight  at  28°C  in  shaker  flasks  containing 
Brain  Heart  Infusion  broth  (Becton,  Dickinson  and  Company  (BD),  East  Rutherford,  NJ),  then 
diluted  directly  from  the  broth  to  the  correct  dosage  based  on  spectrophotometric  optical  den¬ 
sity  readings  (OD600).  The  suspensions  were  prepared  immediately  prior  to  the  administration 
and  kept  on  ice. 

Each  of  the  fasted  and  anesthetized  animals  was  placed  in  dorsal  recumbency  in  the  plethys¬ 
mography  box,  inserting  its  head  through  the  dental  dam  head  port  into  the  head-only  expo¬ 
sure  apparatus  placed  inside  a  Class  3  biosafety  cabinet.  The  bacteria  were  nebulized  using  a 
low-flow  6-jet  Bio  Aerosol  Nebulizing  Generator  (BANG,  CH  Technologies,  Inc.,  Westwood, 
NJ)  with  an  exposure  time  of  between  5  to  15  m.  The  aerosol  particle  stream  was  delivered  at  a 
flow  volume  of  12  1pm  (liter  per  minute)  to  each  anesthetized  monkey  which  was  allowed  to 
breathe  freely.  Real-time  plethysmography  (Skornick  Plethysmography,  CH  Technologies, 

Inc.,  NJ)  was  used  to  measure  tidal,  minute,  and  accumulated  inhaled  volumes  during  exposure 
to  ensure  the  correct  aerosol  dose. 

Aerosolized  bacteria  were  sampled  from  the  head  exposure  chamber  into  an  all-glass  impinger 
containing  20  ml  of  normal  sterile  saline  (Ace  Glass,  Inc.,  Vineland,  NJ)  and  concentrations  were 
confirmed  by  quantitative  bacterial  culture.  Purity  of  the  aerosolized  sample  was  assessed  by  col¬ 
ony  morphology  and  growth  on  Congo  Red-containing  media.  The  target  particle  size,  mean 
mass  aerosol  diameter  (MMAD)  of  between  0.5  and  3  pm,  was  determined  using  a  TSI  Aerosol 
Particle  Sizer  (Model  3321;  TSI,  Inc.,  Shoreview,  MN)  and  by  side-scatter  laser  photometry 
(Microdust  Pro;  Cassella  CEL,  Inc.,  Buffalo,  NY).  Pathogen  dose  was  calculated  using  the  formula 
Dose  =  (C  x  V),  where  C  was  the  concentration  of  viable  pathogen  in  the  exposure  atmosphere, 
and  V  was  the  total  volume  inhaled  based  on  plethysmography  and  exposure  times. 

Blood  collection,  bacteriology  and  biosample  isolation 

From  each  animal,  venous  blood  was  collected  twice  from  the  femoral  vein  (SI  Fig).  During 
each  draw,  samples  were  collected  into  Vacutainer®  CPT  tubes  (BD,  Baltimore,  MD)  and 
PAXgene  tubes  (Qiagen,  Inc.,  Germantown,  MD).  An  aliquot  of  whole  blood  was  plated  on  a 
Congo  red  agar  plate  and  incubated  at  28°C  for  72  h  before  counting  the  colonies.  The  in  vitro 
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Fig  1.  Longitudinal  shift  of  bacteria  load  in  nonhuman  primates’  blood  and  corresponding 
transcriptional  profile.  The  CFU/mL  on  the  Y-axis  (left  side)  is  in  log  scale.  The  circles  indicate  the  bacterial 
load  of  each  animal  and  the  dotted  line  connects  the  mean  bacterial  load  across  the  course  of  infection.  The 
transcriptomic  profile  (Gol -Time)  across  the  course  of  infection  is  also  shown  on  the  right  side  Y-axis.  The 
bars  showing  the  differentially  expressed  transcripts  are  divided  into  up-  and  down-regulated  populations. 


doi:10.1 371/journal,  pone.0151788.g001 


blood  culture  results  are  shown  in  Fig  1  and  S2  Table.  At  the  time  of  euthanasia,  the  body  tem¬ 
perature  of  a  randomly  selected  group  was  recorded.  Since,  the  euthanasia  was  scheduled  till  42 
h  post-exposure,  no  body  temperature  was  recorded  beyond  that  time  point. 

For  pathology,  the  tissues  were  aseptically  collected  and  portions  approximately  1  gm  in 
weight  were  removed.  These  were  placed  into  properly  labeled,  disposable,  sterile,  sample  bags 
that  were  then  heat-sealed,  weighed,  and  placed  on  ice.  Eventually,  these  samples  were  dissoci¬ 
ated  using  a  hand-held  tissue  homogenizer  (Tissue  Tearor,  Inc.)  for  120  seconds,  serially 
diluted  in  sterile  saline,  and  cultured  on  a  Congo  red  agar  plate,  and  incubated  at  28°C  for  72  h 
for  quantitative  measurement  (S2  Table). 

Blood  collected  in  CPT  tube(s)  was  centrifuged  at  1800  xg  (approximately  2800  rpm  on  a 
Sorvall  RT6000  centrifuge)  for  20  m  at  room  temperature.  Without  disturbing  the  whitish  cell 
layer,  the  clear  plasma  from  the  uppermost  layer  was  transferred  to  a  15  ml  tube  and  then 
stored  at  -80°C.  The  plasma  was  used  for  evaluating  the  blood  chemistry. 

For  oligonucleotide  arrays,  whole  blood  was  collected  in  PAXgene  tubes  (QIAGEN,  MD) 
and  total  RNAs  were  isolated  according  to  manufacturer's  instructions,  followed  by  purifica¬ 
tion  using  the  Min  Elute  Cleanup  Kit  (QIAGEN,  MD)  with  the  vendor’s  protocol.  RNA  quality 
and  quantity  were  determined  using  a  Nanodrop  ND-1000  spectrometer  (Thermo  Fisher,  Wil¬ 
mington,  DE).  RNA  from  each  sample  was  assessed  for  purity  by  A26o/28o  ratios.  To  determine 
the  RNA  integrity,  we  analyzed  1  pi  of  eluate  from  each  sample  with  the  Agilent  2100  Bioanaly¬ 
zer  according  to  the  manufacturer’s  instructions. 


Microarray  hybridization  and  post-processing 

Custom  Rhesus  Monkey  Oligo  Microarray  slides  (Catalog  number:  G2519F-015421)  contain¬ 
ing  -40,000  probes  were  obtained  from  Agilent  Technologies  (Agilent  Technologies,  CA).  The 
Agilent  Quick  Amp  Labeling  Kit  (Agilent  Technologies,  CA)  was  used  to  label  BioChain  Mon¬ 
key  universal  reference  RNA  (BioChain  Institute,  Inc.,  Newark,  CA)  with  Cy3  dye  and  tran¬ 
scribed  mRNA  (from  either  control  or  Y.  pestis  exposed  samples)  with  Cy5  dye.  cDNA  was 
hybridized  and  scanned,  and  microarray  images  were  visualized  and  Lowess  normalized  using 
Feature  Extraction  (Agilent  Technologies,  CA);  the  resulting  data  was  analyzed  using  Gene- 
Spring  software  (Agilent  Technologies,  CA). 

The  results  are  available  online  (http://www.ncbi.nlm.nih.gov/geo/),  GEO  ID:  GSE63560. 
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Biological  annotation  and  functional  analysis 

Significantly  altered  transcripts  were  identified  by  two  methods.  Transcripts  altered  at  individ¬ 
ual  time  points  (moderated  t-test  p  <  0.05)  were  pooled  together  and  the  collection  of  tran¬ 
scripts  was  defined  as  Genes  of  Interest-  Time  Course  Differentiation  (Gol  -Time)  (S3  Table  and 
Fig  1  and  S2A  Fig).  This  set  essentially  enlisted  those  transcripts  which  were  significantly  regu¬ 
lated  in  at  least  one  of  the  eight  time  points  investigated  herein. 

Furthermore,  we  segregated  Gol  -Time  into  two  segments,  namely  Early  and  Late  phases  of 
infection  based  on  the  chronological  order  of  blood  draw.  The  Early  phase  of  infection  enlisted 
all  the  transcripts  significantly  altered  in  any  of  the  5  time  points  from  45  m  to  18  h.  Likewise, 
the  Late  phase  of  infection  enlisted  all  the  transcripts  significantly  altered  beyond  18  h  post¬ 
exposure. 

Based  on  the  expression  dynamics,  we  also  identified  two  additional  subsets  from  Gol- 
Time,  namely  Gol -TimeUp  and  Gol-TimeDown  showing  transient  pattern  of  gene  expressions 
correlated  with  the  dynamics  of  pathogenesis.  For  this  purpose,  TimeClust  was  used  for  initial 
filtering  (simplification)  with  Self  Organization  Map  (SOM)  and  Bayesian  Clustering  (S4A 
Fig).  SOM  is  an  unsupervised  neural  network  algorithm  with  prefixed  nodes  (map-units)  dis¬ 
played  on  a  two-dimensional  map  [25].  Since  the  number  of  genes  was  very  large,  we  used  the 
SOM  approach  as  a  preprocessing  step  to  reduce  the  dimensionality  of  the  clustering  problem; 
and  consequently  to  reduce  the  computational  burden  of  the  subsequent  Bayesian  clustering. 
Bayesian  clustering  uses  a  stochastic  population  model  to  represent  the  gene  clusters,  where 
each  cluster  is  assumed  to  be  different  from  the  other  due  to  the  random  effects  (individual  var¬ 
iability).  The  default  setting  of  random  walk  algorithm  scripted  in  TimeClust  program  [26] 
was  used  without  changing  the  value  of  the  highest  Bayesian  Information  Criterion  score  and 
without  using  the  heuristic  strategy  in  order  to  minimize  the  risk  of  a  suboptimal  Bayesian  out¬ 
put  [27]. 

Here,  Gol-TimeLJp  lists  the  transcripts  showing  a  gradual  pattern  of  up-regulation  during 
the  time  course.  In  this  gene  cluster,  individual  gene  expressions  changed  from  down-regula¬ 
tion  during  Early  time  points  to  up-regulation  in  the  Late  phase.  Contrastingly,  Gol-Time¬ 
Down  lists  the  transcripts  showing  a  gradual  trend  of  down-regulation  during  the  same 
temporal  trajectory.  In  this  gene  cluster,  individual  gene  expressions  changed  from  up-regula- 
tion  during  Early  time  points  to  downregulation  in  the  Late  phase.  S4A  Fig  shows  Gol-TimeUp 
and  Gol-TimeDown  hierarchical  clustering  computed  by  Euclidian  distance  algorithm. 

Ingenuity  Pathways  Analysis  (IP A,  QIAGEN,  Inc.,  CA)  identified  those  canonical  pathways, 
which  were  enriched  by  either  up-  or  down-regulated  genomic  pool  (Fig  2  and  S5  Table). 

In  a  parallel  analysis,  we  clustered  all  the  blood  samples  drawn  across  the  time  points  to 
form  two  groups,  namely  the  Y.  pestis  exposed  group  and  the  unexposed  (control)  group.  The 
rationale  behind  this  complementary  analysis  was  to  identify  the  molecular  signatures  associ¬ 
ated  with  Y.  pestis  infection  irrespective  of  the  time-dependent  response.  The  transcripts  signif¬ 
icantly  altered  (moderated  t-test  p  <  0.05)  between  the  exposed  and  control  groups  were 
defined  as  Genes  of  Interest-Infectious  vs.  Control  (Gol -Inf)  (S4  Table).  S2B  Fig  shows  the 
hierarchical  clustering  computed  by  Euclidian  distance  algorithm.  The  Molecular  Activity  Pre¬ 
dictor  (MAP)  tool  of  Ingenuity  Pathway  Analysis  (IPA,  QIAGEN,  Inc.,  CA)  was  used  to  predict 
the  operative  modes  of  the  networks,  i.e.,  the  activated  vs.  inhibited  status  of  the  networks  sig¬ 
nificantly  enriched  by  GoI-/«/(Fig  3  and  S3  Fig). 

Quantitative  Real-time  PCR 

A  panel  of  genes  mined  from  our  microarray  result  was  validated  using  a  real-time  quantitative 
PCR  (qPCR)  assay  as  described  previously  [28].  The  selection  of  these  genes  was  justified  by  a 


PLOS  ONE  |  DOI:1 0.1 371 /journal. pone. 01 51 788  March  22,  201 6 


6/21 


PLOS 


ONE 


Progression  Dynamics  of  Pneumonic  Plague  in  NHP  Model 


□  Transcriptional  Repression 

□  FXR/RXR  Activation 

□  Apoptosis  Signaling 

□  TLR  Signaling 

□  LXR/RXR  Activation 

□  Coagulation  System 

□  Oxidative  Phosphorylation 

□  Flypercytokinemia/  Flyperchemokinemia 

□  T  Cell  Receptor  Signaling 

□  Acute  Phase  Response  Signaling 

| 

b 

< 


□  Acute  Phase  Response  Signaling 

□  T  Cell  and  B  Cell  Signaling 

□  TREM1  Signaling 

□  TREM1  Signaling 

□  NF-kB  ,  TLR  Signaling 
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□  Axonal  Guidance  Signaling 

□  P53,  mTOR  Signaling 

□  Relaxin  Signaling 


□  FXR/RXR  Activation 

□  IL-1  Signaling 

□  GR,  GABA  Receptor  Signaling 
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□  NF-kB,  MAPK,  Gap  Junction  Signaling 

□  Axonal  Guidance  Signaling 


Fig  2.  Potentially  activated  and  inhibited  Gol -Time  enriched  pathways  across  the  time  course.  The  canonical  pathways  enriched  by  up-  or  down- 
regulated  genes  differentially  regulated  at  45  m-6  h,  9  h-1 2  h,  1 8  h-24  h  and  32  h-42  h  post-exposure  are  listed.  The  pathways  constituted  by  the  up-regulated 
genes  are  listed  in  the  “Activated”  section;  the  pathways  constituted  by  the  down-regulated  genes  are  listed  in  the  “Inhibited”  section.  The  functions 
documented  more  than  once  are  in  italics. 

doi:10.1371/journal.pone.0151788.g002 


brief  literature  survey  (Table  1);  those  genes,  which  showed  involvement  with  the  functions 
identified  in  the  present  study  (Results  and  Discussion  of  present  project)  and  past  findings 
(listed  in  the  tables),  were  selected  for  validation.  The  primers  of  the  selected  genes  were  avail¬ 
able  in  catalogued  assays  and  were  purchased  as  a  custom  array  plate  from  SAB  Biosciences 
(Fredrick,  MD).  qPCR  was  carried  out  using  the  RT2  Profiler  PCR  Array  System  (SA  Biosci¬ 
ences,  MD)  on  the  BioRad  CFX  system  (BioRad  Laboratories,  Inc.,  Hercules,  CA)  following 
manufacturers’  instructions.  Two  housekeeping  genes  (ACTB  and  GAPDH)  were  used  for  nor¬ 
malization.  Data  was  analyzed  using  the  Excel-based  RT2  profiler  PCR  Array  Data  Analysis 
using  the  2AAct  method  (Fig  4).  qPCR  outputs  were  compared  to  the  array  results.  We  consid¬ 
ered  these  values  mutually  validated  if  the  outcomes  from  the  two  assay  platforms  showed  reg¬ 
ulation  in  same  direction,  i.e.,  either  both  up-  or  down-regulated  with  -fold  change  greater 
than  ±1.5. 

Blood  Chemistry 

Plasma  samples  obtained  from  the  infected  and  uninfected  nonhuman  primates  (NHP)  were 
assessed  for  chemical  parameters.  These  include:  glucose,  blood  urea  nitrogen  (BUN), 
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CCL3L1/CCL3L3 


Fig  3.  Networks  of  interest.  The  average  fold  changes  of  Early  (45m-1 8h)  and  Late  (24h-42h)  time  groups  are  presented  in  the  right  and  left  boxes, 
respectively.  The  color  scale  presented  in  the  bottom  right  corner  shows  the  gradient  of  fold  change  with  red  upregulated,  blue  downregulated  and  yellow  no 
change  in  gene  expressions.  (A)  Network  of  Stabilization  of  Microtubules.  (B)  Network  profile  of  Ubiquitination. 

doi:  1 0. 1 371  /journal  .pone.0 1 51 788.g003 


creatinine,  sodium,  potassium,  chloride,  carbon  dioxide,  calcium,  phosphorus,  cholesterol,  tri¬ 
glycerides,  total  protein,  albumin,  aspartate  aminotransferase  (Asp  At),  alanine  transaminase 
(ALT),  lactate  dehydrogenase  (LDH),  creatine  kinase  (CK),  alkaline  phosphatase  (ALKP), 
gamma-glutamyl  transpeptidase  (GGT),  and  total  bilirubin.  In  order  to  ensure  safe  transfer  of 
the  samples  out  of  the  BSL-III  facility,  the  plasma  was  first  filtered  through  a  0.2  micron  filter, 
followed  by  screening  of  the  samples  for  bacterial  growth.  Once  all  samples  were  detected  as 
pathogen-free,  300  pi  samples  were  analyzed  (Fig  5).  The  blood  drawn  24  h  before  the  aerosol 
challenge  was  the  baseline  (dotted  horizontal  line,  Fig  5),  which  was  pairwise  (Welch’s  t-test 
p  <  0.05)  contrasted  against  two  post-exposure  clusters,  namely  Early  (45  m-18  h)  and  Late 
(24  h-42  h). 

Results 

Twenty-six  African  green  monkeys  were  subjected  to  aerosol  challenge  with  Y.  pestis  strain 
C092.  Particular  attention  was  paid  to  ensure  equivalent  delivery  of  pathogenic  load.  Blood 
was  drawn  24  h  before  the  exposure  (baseline  control)  and  at  eight  chronologic  time  points  (45 
m,  6  h,  9  h,  12  h,  18  h,  24  h,  32  h  and  42  h)  post-exposure.  Three  control  animals  were  sham- 
exposed  to  saline  under  identical  conditions  similar  to  that  of  Y.  pestis  exposure.  Healthy  rear¬ 
ing  of  those  animals  in  parallel  to  the  exposed  animals  attested  to  the  sanitation  of  their  hous¬ 
ing  environment.  In  addition,  two  infected  controls  were  reared  for  lethality  test  meeting  end 
points  78  h  post-exposure.  The  average  MMAD  ±  geometric  standard  deviation  of  the  gener¬ 
ated  aerosol  particles  was  1.03  pm  ±  1.46,  and  based  on  our  calculation,  the  inhaled  Y.  pestis 
exposure  for  each  animal  was  in  a  range  of  exposures  between  0.33  x  10s  and  3.55  x  106  CFU. 
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Table  1 .  A  brief  literature  survey  underlining  the  relevance  of  the  genomic  pool  selected  for  qPCR. 


Gene  Symbol/  Name 

Description 

IFNG/  Interferon  gamma 

In  the  mouse  pneumonic  plague  model,  IFNG  activation 
was  observed  as  the  consequence  of  pathogenesis  [9]. 

At  the  transcriptomic  level,  IFNG  was  elevated  in  mouse 
models  with  exposure  to  Y.  pestis  strain  of  C092  [11],  but 
suppressed  with  exposure  to  an  avirulent  strain  [49]. 

XIAP/  X-linked  inhibitor  of  apoptosis  protein 

In  vitro  retroviral  infection  mediates  XIAP  to  trigger  the  NF- 
kB  signaling  pathway  as  a  pro-survival  loop  and  to  inhibit 
cell  proliferation  [50]. 

ELP2/  Elongator  Acetyltransferase  Complex 
Subunit  2 

This  Elongator  subunit  is  intricately  associated  with 
chromatin  remodeling  [34]  and  mediates  the  pro-survival 
and  pro-apoptotic  pathways  [51], 

UBE2D1/  Ubiquitin-conjugating  enzyme  E2 

D1 

Yersinia  virulence  factor  interacts  with  UBE2D1  to 
suppress  the  NF-kB  signaling  pathway,  and  thereby 
inhibits  the  downstream  proinflammatory  activities  [52]. 

ADAMTS12/  ADAM  Metallopeptidase  With 
Thrombospondin  Type  1  Motif,  12 

ADAMs  play  key  roles  in  initiating  or  terminating 
immunological  processes;  furthermore,  they  assist  in 
cytokine  biosynthesis  and  direct  systemic  inflammation  or 
barrier  immunity  [53]. 

NCR1/  Natural  cytotoxicity  triggering  receptor 

1 

Initiating  natural  killer  (NK)  cell-mediated  cytolysis,  NCR1 
is  considered  a  typical  marker  for  NK  cells  [54]. 

SOCS1/  Suppressor  of  cytokine  signaling  1 

SOCS1  regulates  the  cytokine  signaling  via  negative 
feedback  loop;  hence,  it  plays  a  critical  role  in  combating 
endotoxic  challenges  [55]. 

NEDD4/  E3  Ubiquitin-Protein  Ligase 

NEDD4  enhances  ubiquitination  and  proteolytic 
degradation  of  many  members  of  the  toll-like  receptor 
channel,  and  thereby  controls  the  early  pathogenic 
recognition  process  [46], 

ALDH1  LI/  Aldehyde  Dehydrogenase  1 

Family,  Member  LI 

Overexpression  of  ALDPI  family  members  helps  to 
ameliorate  oxidative  stress  [56]. 

ALDPI  family  members  control  cellular  motility  [57]. 

IL6/  Interleukin  6 

In  vivo  exposure  of  LPS  surges  the  IL-6  load  in  equine 
whole  blood  [58]. 

A  study  of  the  patients  in  surgical  units  found  IL-6  as  an 
early  tool  to  diagnose  severe  sepsis  [59], 

doi:10.1371/journal.pone.0151788.t001 


The  quantitative  measurements  of  Y.  pestis  present  in  the  exposed  animals’  blood  and  urine 
samples  and  nine  tissues,  including  axillary  lymph  node,  submandibular  lymph  node,  inguinal 
lymph  node,  mediastinal  lymph  node,  mesenteric  lymph  node,  spleen,  liver,  lung  and  kidney, 
were  evaluated  (Fig  1  and  S2  Table)  based  on  the  geometric  mean  of  CFU/ml  of  the  dilution 
factor  and  colony  counts.  Bacterial  colonies  were  detectable  at  9  h  post-exposure  among  a  sub¬ 
set  of  the  infected  cohort.  Some  of  the  results  reported  herein  failed  to  pass  the  threshold  value 
(>  200  CFU/ml);  nevertheless,  they  are  presented,  since  we  noticed  viable  colonies  of  the  inoc¬ 
ulated  bacteria.  Notably,  we  found  no  correlation  of  the  bacterial  enumeration  in  blood  with 
the  inhaled  dose  for  individual  NHPs.  Bacteria  were  detected  in  blood  9  h  after  exposure  in  1  of 
3  animals,  2  of  4  at  12  h,  1  of  3  at  18  h,  3  of  4  at  24  h,  3  of  3  at  32  h,  and  3  of  3  at  42  h.  Even  so, 
liver  and  lungs  showed  very  early  and  consistent  load  of  infection  from  6  h  post-exposure. 

Regarding  other  pathological  data,  we  noticed  minimum  changes  in  the  body  temperature 
between  6  h  and  42  h  post-exposure  with  a  range  of  34.2°C  to  37.3°C  and  the  average  tempera¬ 
ture  was  36.01±0.8°C.  Much  of  the  temperature  differences  were  attributed  to  the  individual 
variances;  maximum  temperature  change  in  individual  NHP  recorded  in  this  study  was  1.3°C 
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Fig  4.  qPCR  results  of  the  selected  genes  arranged  in  the  descending  order  of  validation  (%  in 
parenthesis  beside  the  gene  name  shows  the  number  of  genes  -fold  changed  in  same  direction,  >±1 .5) 
between  Array  (A,  top  row)  vs.  qPCR  (P,  bottom  row)  across  the  time  course,  two  from  Early  (6  h  and 
1 2  h)  and  two  from  Late  time  groups  (32  h  and  42  h).  The  color  scale  presented  in  the  bottom  right  corner 
shows  the  gradient  of  fold  change  with  red  upregulated,  blue  downregulated  and  yellow  no  change  in  gene 
expressions. 


doi:10.1 371/journal,  pone.  0151788.g004 


with  average  change  in  individual  NHP  was  only  0.1 8°C.  Those  are  in  agreement  with  an  ear¬ 
lier  report  [12]. 


Transcriptomic  analysis 

In  comparison  with  the  baseline  defined  by  the  24  h  pre-exposed  blood  samples,  the  differen¬ 
tially  regulated  transcripts  at  the  post-exposure  time  points  were  identified  using  moderated  t- 
test  ( p  <  0.05).  Fig  1  and  S2A  Fig  show  that  the  number  of  significantly  differentially  expressed 
transcripts  gradually  increased  from  45  m  post-exposure,  reaching  a  maximum  during  18  h 
and  24  h,  with  a  subsequent  decrease.  Together,  there  were  10,059  unique  transcripts  enlisted 
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CK  LDH  AspAt  BUN 


Fig  5.  Longitudinal  profile  of  plasma  chemistry.  The  expressions  of  selected  chemicals  (Exposed/ 
Controls)  are  shown.  The  values  measured  from  45  m  to  1 8  h  were  pooled  together  as  the  Early  group  (white 
bar)  and  those  from  24  h  to  42  h  were  clustered  as  the  Late  group  (shaded  bar).  The  horizontal  broken  line 
depicts  the  control  value  calculated  from  the  blood  samples  drawn  24  h  before  aerosol  exposure.  Welch’s  t- 
test  computed:  #  0.1  <  p  <  0.05,  *  p  <  0.05,  **  p  <  0.001 ,  ***  p  <  0.0001 . 

doi:10.1 371/journal,  pone.0151788.g005 


by  Gol -Time,  but  no  transcripts  were  found  conserved  across  all  the  time  points  investigated 
herein. 

Furthermore,  Gol- Time  (S2A  Fig  and  S3  Table)  was  segregated  in  chronologic  order.  As 
explained  earlier,  Early  phase  of  infection  covered  from  45  m  to  18  h  post-exposure  time 
points,  while  the  Late  phase  of  infection  covered  the  remainder,  starting  from  24  h  post-expo- 
sure.  Making  the  demarcation  between  18  h  and  24  h  was  justified  by  the  bell-shaped  transcrip¬ 
tional  expression  profile  (Fig  1  and  S2A  Fig),  which  peaked  between  18  h  and  24  h.  There  were 
6,495  and  6,020  transcripts  altered  during  the  Early  and  Late  phases  of  infection,  respectively, 
while  2,456  transcripts  were  conserved  between  these  two  phases  of  infection.  Using  the  SOM 
algorithm,  we  further  sorted  Gol- Time  to  mine  the  transcripts  showing  consistent  temporal 
changes  either  towards  up-regulation  (138  transcripts)  or  down -regulation  (225  transcripts) 

(S4  Fig)  designated  as  Gol-TimeUp  and  Gol -TimeDown,  respectively. 

In  a  parallel  analysis,  we  identified  1,274  transcripts  defined  as  GoI-/«/(S2B  Fig  and  S4 
Table)  that  were  differentially  regulated  between  the  exposed  and  baseline  values;  the  latter  was 
derived  from  the  blood  drawn  24  h  pre-exposure.  GoI-/n/was  essentially  the  collection  of  all 
the  transcripts  identified  in  the  blood  samples  drawn  across  the  entire  time  course  post  Y.  pestis 
exposure.  Interestingly,  a  considerable  number  of  Gol-Inf  transcripts  (82%)  were  found  con¬ 
served  with  Gol- Time  (S2C  Fig). 


Functional  analysis  of  Gol -Time 

A  brief  overview  of  the  canonical  functional  analysis  of  Gol-Time  is  depicted  in  Fig  2  and  the 
details  are  listed  in  S5  Table.  In  addition,  parallel  functional  analyses  of  Gol-TimeUp  and  Gol- 
TimeDown  were  performed  to  elucidate  time-dependent  modifications  of  the  functional  activi¬ 
ties.  S4B  Fig  depicts  the  hierarchical  clustering  of  25  transcripts  from  Gol-TimeUp  associated 
with  the  most  enriched  functional  modules.  These  genes  were  further  clustered  into  three  sub¬ 
groups  (A up,  B up  and  Cup)  based  on  their  biological  interactions  mapped  from  literature  min¬ 
ing.  S4B  Fig  also  displays  the  hierarchical  clustering  of  86  transcripts  from  Gol-TimeDown 
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associated  with  the  most  enriched  modular  interaction  networks.  The  corresponding  networks 
are  shown  in  S4C,  S4D  Fig. 

Functional  analysis  of  Go\-lnf 

Functional  annotation  of  Go  I -/n/ (S3  Fig)  demonstrated  a  significant  enrichment  in  apoptosis 
(S5A  Fig),  stabilization  of  microtubules  (Fig  3A),  recruitment  of  inflammatory  leukocytes  (S5B 
Fig),  formation  of  cellular  inclusion  bodies,  synthesis  of  extracellular  matrix,  and  myocardial 
and  liver  failure.  Protein  ubiquitination  signaling  (Fig  3B),  G2/M  DNA  damage  checkpoint, 
growth  hormone  signaling,  and  p53  signaling  pathways  were  the  top  enriched  canonical 
pathways. 

The  activation  or  inhibition  dynamics  of  the  networks  across  the  infection  period  were  pre¬ 
dicted  by  MAP  (IP  A,  CA)  with  the  relationships  scored  in  silico  based  on  the  regulations  of  the 
molecules  connected  to  the  relevant  annotation  (S3  Fig).  MAP  predicted  a  consistent  activation 
of  apoptosis  encompassing  47  transcripts.  Fourteen  and  eight  apoptotic  transcripts  were  associ¬ 
ated  with  DNA  repair  and  chromatin  condensation,  respectively,  although  with  a  limited  pre¬ 
dictive  power  due  to  small  sample  size.  MAP  further  predicated  temporally  consistent 
activation  of  the  recruitment  of  inflammatory  leukocytes  and  ubiquitination,  and  potential 
inhibition  of  the  microtubule  stabilization  across  the  time  course. 

Validation  using  qPCR  results  (Fig  4) 

The  high  throughput  microarray  data  was  validated  by  real-time  qPCR.  A  brief  literature  sur¬ 
vey  outlining  the  relevance  of  the  ten  genes  selected  for  this  study  is  described  in  Table  1.  For 
the  qPCR  study,  two  time  points  were  selected  from  Early  (6  h  and  12  h)  and  Late  (32  h  and  42 
h)  phases  of  infection,  respectively.  Fig  4  depicts  the  comparative  analysis  of  the  oligonucleo¬ 
tide  microarray  result  and  the  qPCR  result.  The  number  of  transcripts  showing  similar  regula¬ 
tions  in  both  platforms  was  calculated.  Six  transcripts  showed  similar  expression  in  three  out  of 
four  time  points,  and  four  transcripts  showed  similar  -fold  changes  in  two  out  of  four  time 
points. 

Blood  Chemistry  (Fig  5) 

The  blood-based  evaluation  was  limited  to  the  blood  chemistry  of  the  serum.  A  detailed  proto¬ 
col  was  described  that  highlights  our  effort  to  maintain  the  serum  pathogen-free.  For  safety  rea¬ 
sons,  we  could  not  inspect  the  basic  cellularity  of  blood  potentially  carrying  BSL-3  pathogens. 

It  was  a  drawback  of  this  study. 

To  compute  the  blood  chemistry,  we  clustered  the  exposed  animals  into  two  groups  in  an 
attempt  to  enhance  the  statistical  power  of  the  analysis.  As  noted  previously,  the  Early  phase  of 
infection  consisted  of  the  samples  collected  over  the  time  course  from  45  m  to  1 8  h,  whereas 
the  Late  phase  of  infection  consisted  of  the  samples  collected  from  24  h  to  42  h.  Comparisons 
of  the  Early  or  Late  phase  of  infection  versus  the  baseline  controls  (pre-exposure  blood) 
showed  CK  and  AspAt  consistently  elevated.  BUN  and  LDH  were  significantly  elevated  in  the 
Early  phase  of  infection,  but  these  returned  to  the  baseline  during  the  Late  phase  of  infection. 

Discussion 

We  studied  the  longitudinal  dynamics  of  perturbation  of  the  blood  transcriptome  of  the  Afri¬ 
can  green  monkey  ( Chlorocebus  aethiops)  as  a  consequence  of  aerosol  exposure  to  Y.  pestis. 

The  baseline  of  the  study  was  defined  by  the  blood  samples  drawn  24  h  before  exposing  the 
monkeys  to  Y.  pestis;  the  post-exposure  blood  samples  were  drawn  at  eight  sequential  time 
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points  (SI  Fig).  Two  of  the  treated  animals  were  allowed  to  live  past  42  h  and  met  end  points  at 
nearly  78  h  post-exposure. 

The  lethality  test  showed  that  the  pathogen  load  used  in  the  present  study,  if  left  untreated, 
can  cause  death  within  78  h  post- exposure.  Therefore,  to  characterize  the  early  dynamics  of 
pathogenesis,  we  investigated  the  transcriptomic  profile  till  the  half-way  through  the  survival 
curve,  which  was  approximately  42  h  post-exposure. 

The  animals  showed  no  elevation  of  body  temperature  till  42  h  post  infection.  No  viable  bac¬ 
terial  presence  was  recorded  until  9  h  post-exposure;  in  fact,  only  a  fraction  of  infected  NHPs 
manifested  bacterial  colonization  in  the  blood  prior  to  24  h  post-exposure,  which  was  in  agree¬ 
ment  with  the  previous  report  [12].  Further  investigation  found  bacterial  loads  in  kidney,  lungs 
and  liver  at  much  earlier  time  points  (S2  Table).  It  could  be  argued  that  the  pathogens  were 
possibly  nested  inside  and  transported  via  the  macrophages  and  other  blood  cells  during  this 
Early  phase  of  infection,  thereby  evading  in  vitro  detection  in  blood  [11].  Other  possibilities 
such  as  potential  delayed  emergence  of  bacteria  from  the  pulmonary  tract  and/or  the  lymph 
nodes  could  not  be  overruled,  and  this  underlines  the  necessity  of  further  investigation. 

The  molecular  perturbations  were  noted  at  45  m  post-exposure  (Fig  1),  much  earlier  than 
the  emergence  of  the  bacterial  colonies  in  blood;  although  we  cannot  overrule  the  potential 
contributions  of  certain  external  factors,  such  as  anesthesia,  in  changing  the  early  genomic  pro¬ 
file.  In  addition,  this  study  presented  the  transcriptomic  regulation  averaged  over  all  the  blood 
cell  types.  Cell- specific  data  could  produce  a  wealth  of  information,  but  this  was  beyond  the 
scope  of  the  present  approach. 

Oligonucleotide  microarray  data  was  analyzed  in  two  complementary  fashions.  First,  we 
probed  the  transcriptomic  profile  at  individual  post-exposure  time  points  (Gol-Time).  Second, 
we  contrasted  the  entire  infected  cohort,  irrespective  of  their  infection  duration,  against  their 
baseline  derived  from  the  blood  drawn  24  h  pre-exposure  (GOI -Inf).  It  is  important  to  note 
that  the  present  result  is  limited  by  probing  the  Rhesus  Macaques  (Macaca  mulatta)  cDNA 
array,  but  rhesus  macaque  was  phylogenetically  diverged  from  African  Green  Monkey  10-15 
million  years  ago.  This  phylogenetic  divergence  could  bias  the  results  regarding  the  most  rap- 
idly-evolving  immune  gene  families,  such  as  variable  cell-surface  ligands  or  variable  receptors; 
some  of  those  have  been  identified  in  the  present  study.  To  mitigate  the  potential  bias  due  to 
the  evolutionary  divergence,  a  next  generation  sequencing  study  is  recommended. 

Characterizing  GOI -Time,  we  observed  that  the  number  of  transcripts  gradually  increased 
until  24  h,  followed  by  a  steady  reduction.  The  initial  rise  of  transcriptomic  numbers  was 
potentially  due  to  the  accumulative  influences  of  the  pathogenic  infection.  However,  at  later 
phases  of  infection,  simultaneous  activation  of  many  intertwined  networks  potentially  trig¬ 
gered  conflicting  feedbacks  to  the  transcriptional  activities.  As  a  result,  the  number  of  differen¬ 
tially  expressed  genes  was  reduced  among  the  moribund  animals.  Similar  temporal 
distribution  of  the  genomic  profile  depicting  a  bell-shaped  curve  was  observed  previously  in  a 
study  of  Mycobacterium  tuberculosis  in  vivo  [29] . 

Gol-Time  elucidated  valuable  information  about  the  molecular  shift  as  the  function  of  the 
time  lapse  after  the  pathogenic  challenge.  This  strategy  was,  however,  particularly  vulnerable  to 
the  limited  statistical  power.  The  sample  size  of  each  time  point  was  limited  to  three  or  four. 
Therefore,  we  complemented  the  analysis  with  GoI-/w/ as  the  genomic  profile  of  all  post- 
exposed  26  NHPs  were  clustered,  irrespective  of  their  infection  durations,  and  contrasted 
against  the  baseline  derived  from  the  blood  drawn  24  h  pre-exposure.  With  this  larger  sample 
size,  we  expected  a  better  statistical  confidence  in  elucidating  the  pathogenic  implications  over 
the  entire  infection  period  while  at  a  cost  of  the  finer  detail  of  time-specific  information.  Inter¬ 
estingly,  around  83%  of  the  transcriptomic  members  of  Gol-Inf  were  also  conserved  in  Gol- 
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Time,  and  they  were  almost  equally  distributed  between  the  two  Early  (45  m  -  18  h)  and  Late 
(24  h — 42  h)  phases  of  infection. 

A  systematic  modulation  of  the  host  immunity  initiated  at  the  onset  of  infection  was  evident 
from  transcriptomic  results.  Supporting  a  previous  study  [8],  Y.pestis  first  suppressed  the 
genes  enriching  networks  associated  with  the  natural  killer  (NK)  cell  signaling  and  phagocyto¬ 
sis  in  macrophages  to  overwhelm  the  innate  immunity  (Fig  2).  A  pool  of  pro-inflammatory 
transcripts,  such  as  SOCSS1,  UBE2D1,  and  ADAMTS12,  was  differentially  regulated  as  shown 
by  qPCR  results  (Fig  4).  Many  cytokines,  such  as  IL-1A,  -11,-18,  and  -27,  and  chemokines, 
such  as  CCL-1,  -14,  -17,  and  -26,  and  CXCL16  showed  gradual  transcriptomic  suppression. 
Early  onset  of  apoptosis  manifested  by  the  early  perturbation  of  JAK/Stat,  PI3K/AKT  and 
MAPK  pathways  (Fig  2)  potentially  played  a  significant  role  in  facilitating  Y.  pestis  evasion  of 
early  pathogenic  recognition  [30,31].  In  fact,  the  active  toll-like  receptor  (TLRs)  signaling  net¬ 
works  (Fig  2)  possibly  failed  to  contain  Y.  pestis.  Further  studies  are  necessary  to  confirm  this 
hypothesis. 

Delayed  suppression  of  pathways  linking  MAPK  and  NF/cB  signaling  potentially  triggered 
rapid  apoptosis  or  necrosis  at  the  later  time  points  (S3  Table)  [32].  In  addition,  mTOR,  AMPK 
and  p53  signaling  networks  were  also  suppressed  beyond  24  h  post-exposure  (S3  Table).  Those 
potentially  facilitated  the  pathogen-promoted  annihilation  of  the  host  defense  [33].  As  a  conse¬ 
quence,  two  biofunctions  closely  linked  with  apoptosis,  namely  the  repair  of  DNA  and  chroma¬ 
tin  condensation,  were  suppressed.  The  qPCR  result  reported  the  consistent  up-regulation  of 
ELP2,  a  mediator  of  the  chromatin  remodeling  process  [34]. 

In  a  parallel  chain  of  events  in  the  Early  phase  of  infection,  activation  of  many  inflammatory 
pathways,  such  as  FXR/RXR,  TREM,  iNOS  and  LXR/RXR  triggered  the  cytokine  production  in 
macrophages  and  oxidative  phosphorylation  (Fig  2).  The  increasing  expression  of  TLR  mole¬ 
cules,  though  failing  in  pathogenic  recognition,  potentially  activated  NFxB  and  thereby  initi¬ 
ated  the  pro-survival  loop  [35].  Subsequent  elevation  of  the  genes  encoding  the  inflammatory 
chemokines  such  as  CCL16,  CCL25  and  CXCL6,  and  cytokines  such  as  IL-2,  IL-8  and  IFNG 
(Fig  4)  in  concert  with  delayed  activation  of  the  coagulation  system  (Fig  2)  are  the  typical  signa¬ 
tures  of  sepsis  [36]  [37]. 

With  the  first  indication  of  an  uncontrolled  surge  of  cytokines  and  chemokines  (activated 
networks  of  hypercytokinemia  and  hyperchemokinemia,  as  shown  in  Fig  2),  24  h  seemed  to  be 
the  time  point  when  Y.  pestis  began  to  overwhelm  the  host  defense  mechanisms.  Beyond  the  24 
h  post-exposure,  Y.  pestis  began  to  perturb  many  cancer-related  networks  affecting  organs 
such  as  pancreas,  breast,  lungs,  liver  and  bone  marrow  (S5  Table).  An  uncontrolled  surge  of 
cytokines  and  chemokines  was  suggested  beyond  24  h  post-infection,  which,  coupled  with 
multi-organ  failures,  could  bring  about  fatality  in  rapid  succession  [38] . 

Supportive  of  past  evidence  [39],  the  present  report  also  indicates  that  Y.  pestis  infection 
potentially  delayed  the  activation  of  adaptive  immunity.  The  networks  associated  with  T-cell 
and  B  cell-activation  were  triggered  beyond  the  24  h  post-exposure  (Fig  2).  These  late 
responses  were  likely  to  be  ineffective  in  containing  the  already  well-proliferated  pathogenesis. 
The  present  study,  in  concurrence  with  many  previous  observations  [40-42],  suggests  value  in 
including  priming  of  cellular  immunity  in  the  next  generation  vaccination  program. 

The  longitudinal  dynamics  of  the  transcriptomic  shift  (Fig  1)  suggestively  contradicts  the 
biphasic  Y.  pestis  progression  pattern  reported  in  the  mouse  model  [9,22],  which  portrayed  an 
early  phase  of  subdued  recruitment  of  inflammatory  markers  turning  into  the  uncontrolled 
discharge  of  cytokines  and  chemokines.  However,  a  recent  study  showed  an  early  response  of  a 
very  different  quality  of  cytokines  in  a  potentially  plague -resistant  strain  of  mouse  [19]. 
Extending  the  knowledge,  the  present  NHP  model  demonstrated  an  Early  pro-inflammatory 
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phase  evidently  mediated  by  active  ubiquitin  and  lymphocytes,  which  eventually  turned  into 
hypercytokinemia  and  hyperchemokinemia. 

In  the  context  of  early  immune  response,  we  draw  attention  to  our  observation  of  hemato¬ 
logical  chemistry  (Fig  5).  The  early  spike  of  BUN  and  LDH  load  potentially  reflects  the  cyto¬ 
kine-mediated  muscular  catabolism  [43]  and  cytotoxicity  [44]  ensuing  immediately  after  the  Y. 
pestis  exposure.  The  early  elevation  of  CK  and  Asp  At  possibly  gives  an  estimation  of  the  scale 
of  rapidity  with  which  Y.  pestis  was  disseminated  to  liver  and  muscle  tissues.  These  observa¬ 
tions  along  with  the  bacterial  load  enumeration  matrix  (S2  Table)  potentially  support  the 
hypothesis  that  the  pathogens  were  transported  via  macrophages  and  other  blood  cells  during 
the  Early  phase  of  infection,  and  thereby  escaped  early  detection  (Fig  1)  [11].  Comprehensive 
transcriptomic  investigation  focused  on  the  tissues  of  interest  is,  therefore,  crucial  in  order  to 
understand  the  route  of  pathogenesis  and  the  pertinent  mechanisms. 

The  role  of  microtubules  in  Y.  pestis  uptake  via  epithelial  cells  was  suggested  previously 
[45].  Epithelial  cells  served  as  an  essential  port  for  the  inhalation  exposure  in  the  present  study. 
Therefore,  the  genes  encoding  the  early  destabilization  of  the  microtubules  (Fig  3A)  could  have 
therapeutic  implications.  Furthermore,  longitudinal  annotations  of  the  transcriptomic  profile 
predicted  a  consistent  elevation  of  ubiquitination  (Fig  3B).  qPCR  results  showed  early  overex¬ 
pression  of  NEDD4,  a  genomic  activator  of  ubiquitination  [46].  The  contribution  of  ubiquiti¬ 
nation  to  cellular  immunity  [47],  particularly  in  mediating  NFkB  activation  [48],  could  be 
explored  further  in  the  quest  for  an  effective  therapy  against  Y.  pestis. 


Conclusions 

The  transcriptomic  analyses  suggest  a  very  rapid  molecular  perturbation  in  blood  was  caused 
by  Y.  pestis  assault  in  African  green  monkeys.  The  host  defense  was  essentially  exhausted,  and 
the  process  of  multi- organ  failure  was  initiated  only  one  day  after  the  aerosol  exposure,  while 
the  body  temperatures  showed  no  elevation.  Therefore,  we  speculate  that  the  molecular  signa¬ 
tures  mined  prior  to  the  24  h  post-exposure  could  have  early  therapeutic  potential. 

During  the  early  phase  of  pathogenesis,  systematic  modulation  of  both  adaptive  and  innate 
immunity  programs  was  mediated  in  parallel  with  activation  of  toll-like  receptor-based  patho¬ 
genic  detection.  Despite  that,  pathogenesis  progressed  unabated.  We  hypothesize  that  the  early 
onset  of  apoptosis  facilitated  the  disease  proliferation  [30,31].  This  hypothesis  is  subject  to  veri¬ 
fication  through  an  independent  study. 

We  also  observed  early  inhibition  of  natural-killer  cell  signaling  networks  causally 
compromising  the  innate  immunity,  in  agreement  with  previous  observations  [39].  The 
delayed  activation  of  adaptive  immunity  would  likely  to  come  out  ineffective  in  combating 
highly-proliferating  pathogens.  Coordinated  priming  of  innate  and  adaptive  immunity  is 
thereby  suggested  as  a  potentially  effective  method  for  next  generation  vaccination  programs 
[40-42].  In  the  context  of  immunological  health,  the  involvement  of  ubiquitin  in  Y.  pestis  path¬ 
ogenesis  is  worth  exploring  in  the  future. 

African  green  monkeys  showed  an  early  surge  of  pro-inflammatory  immune  response  via 
activation  of  lymphocytes,  which  was  a  possible  sign  of  deviation  from  the  biphasic  Y.  pestis 
pathogenesis  model  typically  reported  in  rodents  [9,22].  Since  the  present  NHP  model  is  phylo- 
genetically  closer  to  human,  our  observations  could  be  very  important  in  defining  next  genera¬ 
tion  care;  however,  further  validation  is  necessary. 

Along  this  line,  we  also  must  admit  that  we  found  a  number  of  similarities  between  our 
NHP  model  and  previously  published  murine  models.  Many  of  the  genes  identified  in  the  pres¬ 
ent  study,  such  as  IL2RG,  NdRG,  IFNG,  and  CXCL6,  were  also  regulated  in  the  mouse  pneu¬ 
monic  plague  models  [15,49],  Increased  IFNG  and  IL-6  loads  were  also  reported  in  Y.  pestis- 
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infected  rodent  lungs  [9].  Furthermore,  many  networks  identified  in  the  present  study  are  in 
agreement  with  the  murine  tissue  model  [11],  namely  the  networks  linked  to  the  antigen  pre¬ 
sentation,  cell  death,  cell-mediated  immunity,  inflammatory  response,  hematopoiesis,  acute 
phase  response,  glucocorticoid  receptor  signaling,  DNA  repair,  free  radical  scavenging  and 
LXR/RXR  activation. 

In  summary,  the  present  study  recognizes  a  potential  pattern  of  Y.  pestis  pathogenesis  in 
evasion  of  the  early  host  defense  mechanism.  The  role  of  microtubules,  ubiquitination,  and  the 
contribution  of  early  apoptotic  aggravation  suggested  in  this  context  are  subject  to  further  vali¬ 
dation.  Follow-up  comprehensive  investigations  focused  on  organs  of  interest  are  underway  to 
further  elucidate  the  pathogenesis  and  disease  proliferation. 

Supporting  Information 

51  Fig.  Flow  chart  of  the  study  design.  The  study  began  with  initial  blood  draws  conducted  at 
T0-24  h.  At  T0,  nonhuman  primates  were  exposed  to  Y.  pestis  (marked  as  EXP+).  Final  blood 
draws  were  carried  out  at  eight  different  time  points  marked  as  T0  +  45  m,  6  h,  9  h,  12  h,  18  h, 
24  h,  32  h  and  42  h.  The  number  of  primates  ( N )  used  for  each  time  point  is  indicated  as  well. 
All  of  the  animals  except  the  45  m  cohort  were  euthanized  to  collect  organs  stipulated  for 
another  study.  The  same  cohort  of  animals  was  bled  twice,  at  45  m  and  6  h  post-stress. 

(TIF) 

52  Fig.  Comparative  analysis  of  Gol-Time  and  Gol-Inf.  The  transcriptomic  profile  of  Gol- 
Time:  across  the  time  course.  The  volcano  graph  is  a  plot  of  p-value  [-loglO(p-value)]  against 
-fold  change  [log2(fold  change)].  (B)  Hierarchical  clustering  of  Gol-Inf  across  the  time  course. 
The  color  scale  presented  in  the  bottom  right  corner  shows  the  gradient  of  fold  change  with 
red  upregulated,  blue  downregulated  and  yellow  no  change  in  gene  expressions.  (C)  Relation¬ 
ship  between  Gol-Inf  and  Gol-Time  transcripts.  The  Venn  diagram  shows  the  overlap  of  tran¬ 
script  populations  between  Gol-Inf  and  the  two  sets  of  Gol-Time.  These  two  subsets  of 
Gol-Time  are  the  following:  (i)  the  Early  group  is  a  pool  of  the  transcripts  significantly  altered 
at  45  m,  6  h,  9  h,  12  h  and  18  h;  (ii)  the  Late  group  is  a  pool  of  the  transcripts  significantly 
altered  at  24  h,  32  h  and  42  h. 

(TIF) 

53  Fig.  Longitudinal  dynamics  of  functions  enriched  by  Gol-Inf.  The  functional  regula¬ 
tions  were  computed  by  Molecular  Predicative  Analysis  (MAP),  a  z-score-based  algorithm 
provided  by  IPA.  The  up-  and  down-wards  arrowheads  depict  the  activated  and  inhibited  lev¬ 
els  of  these  functions.  The  non-active  levels  are  shown  by  ‘X’.  The  functions  of  interests  are  (i) 
Apoptosis,  (ii)  Stabilization  of  microtubules,  (iii)  Recruitment  of  inflammatory  leukocytes,  (iv) 
Ubiquitination,  (v)  DNA  repair,  and  (vi)  Chromatin  condensation. 

(TIF) 

54  Fig.  Characterization  of  GOI -TimeUp  and  Gol -TimeDown.  (A)  Hierarchical  clusters  of 
two  subsets  of  Gol-Time  showing  a  significant  gradual  trend  of  shifting  towards  (i)  overexpres¬ 
sion  (Gol -TimeUp)  and  (ii)  suppression  (Gol -TimeDown).  The  gene  sets  were  named  Gol- 
TimeUp  and  Gol -TimeDown,  respectively.  The  color  scale  presented  in  the  bottom  right  corner 
shows  the  gradient  of  fold  change  with  red  upregulated,  blue  downregulated  and  yellow  no 
change  in  gene  expressions.  (B)  Figure  S4(B).  Characterization  of  GO  I  -  Time  Up  and  Gol-Time- 
Down.  Hierarchical  clustering  of  the  focus  list  from  Gol-TimeUp  and  Gol  -TimeDown-  The  left 
hierarchical  cluster  depicts  the  subsets  from  Gol -TimeUp  showing  the  most  enriched  modular 
interactions.  They  are  grouped  in  three  sub-clusters  (Aup-Cup)  and  the  related  functions  are 
listed  (in  bullets)  to  its  right.  During  the  infection  lifecycle,  these  functions  were  consistently 
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enriched  by  gradually  overexpressing  genes.  The  right  hierarchical  cluster  depicts  the  subsets 
from  Gol-TimeDown  showing  the  most  enriched  modular  interactions.  They  are  grouped  in 
fives  (Ajown  -EdoWn)  and  the  related  functions  are  listed  (in  bullets)  to  its  right.  During  the 
infection  lifecycle,  these  functions  were  consistently  enriched  by  gradually  suppressing  genes. 
The  color  scale  presented  in  the  bottom  right  corner  shows  the  gradient  of  fold  change  red 
upregulated,  blue  downregulated  and  yellow  no  change  in  gene  expressions.  (C)  Network  of 
the  focus  subset  from  Gol -TimeUp.  As  shown  in  Figure  S4B,  the  genes  are  clustered  into  the 
group  identified  as  A-Cup.(D)  Network  of  the  focus  subset  from  Gol-TimeDown.  As  shown  in 
Figure  S4B,  the  genes  are  clustered  into  the  groups  identified  as  A-E<£0wn. 

(TIF) 

S5  Fig.  Networks  of  interest.  The  candidate  genes  are  colored  by  their  regulation  at  24  h.  The 
color  scale  presented  in  the  bottom  right  corner  shows  the  gradient  of  fold  change  with  red 
upregulated,  blue  downregulated  and  yellow  no  change  in  gene  expressions.  (A).  Network  of 
the  apoptosis  profile.  (B).  Network  of  the  recruitment  of  inflammatory  leukocytes. 

(TIF) 

51  Table.  The  lethality  endpoints  approved  under  IACUC  approved  protocol  P001-08:  Sys¬ 
tems  biology  studies  to  identify  host  indicators  /therapeutic  targets  at  early  time  periods 
post  exposure  to  Y.  pestis  (C092)  in  African  Green  Monkeys  (Chlorocebus  aethiops).  Any 

animal  found  to  present  with  3  or  more  of  the  clinical  signs  or  to  be  moribund,  cachectic,  or 
unable  to  obtain  food  or  water  will  be  judged  to  be  at  or  near  the  endpoint  and  humanely 
euthanized 
(XLSX) 

52  Table.  Bacterial  load  enumerated  from  different  organs  of  infected  animals.  *Not  appli¬ 
cable;  **No  animals  were  euthanized  at  45m;  therefore  no  bacterial  loads  were  detected  in  tis¬ 
sues.  These  same  animals  were  euthanized  at  the  6h  time  point.  ***One  animal  died  at  72h  and 
the  other  one  was  euthanized  at  78h;  so  all  the  values  in  bacterial  loads  were  enumerated  at  78h 
euthanasia  time  points. 

(XLSX) 

53  Table.  The  Gol -Time  regulation  profile.  This  is  the  list  of  genes  significantly  regulated 
at  individual  time  points.  The  values  are  in  log2(fold  change). 

(XLSX) 

54  Table.  The  Gol-Inf  regulation  profile.  This  is  the  list  of  genes  identified  by  contrasting  all 
the  post-exposure  samples  (irrespective  of  time  of  infection)  vs.  the  respective  controls  (bloods 
drawn  before  the  aerosol  challenge).  The  values  are  in  log2(fold  change). 

(XLSX) 

55  Table.  The  canonical  pathways  enriched  by  Gol -Time. 

(XLSX) 

56  Table.  Genes  associated  with  the  functions  significantly  enriched  by  Gol-Inf.  The  values 
are  in  log2(fold  change). 

(XLSX) 
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